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FOREWORD 

This  document  replaces  NEL  Technical  Memorandum  600 
Kectetc  o,  Power  Speer™  md  Cross  Analysis  by  Digital  Method s,  15 
■  pul  I9W.  which  was  a  compilation  of  lecture  notes  on  power-spectrum 
analvsts  by  Dr.  E.  E.  Gossard.  The  information  has  been  reorganized 

added'1"""1111  lll"St,aUo"s  a,ld  Partinent  literature  references  have  been 
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INTRODUCTION 


Complicated  interactions  of  physical  processes  often  occur  in 
nature  and  their  recorded  measurements  appear  random  or  pseudorandom 
.in  character.  Such  records  occur  especially  often  in  dealing  with  atmos¬ 
pheric  turbulence  or  broadband  wave  phenomena.  The  medium  is  usually 
far  too  complicated  to  be  predicted  or  even  described  in  detail. 

The  autocorrelation  function  was  first  used  as  an  independent 
variable  in  geophysical  problems  by  G.  I.  Taylor1  in  his  studies  of  turbu¬ 
lence.  It  represented  an  important  advance  in  studying  complicated 
physical  processes  and  forms  the  basis  of  modern  turbulence  theory.  The 
notions  of  Taylor  were  extended  and  explored  by  Von  Karmen  and  others. 
The  mathematical  theory  relating  correlation  functions  and  spectra  has 
been  examined  in  detail  by  Wiener.2 

The  usefulness  of  the  cross  or  autocorrelation  function  is  based 
on  the  fact  that  the  statistical  character  of  an  extremely  complicated  re¬ 
cord  may  often  be  described  by  a  ratter  simple  correlation  function 
(Gaussian,  Markov,  etc.).  While  details  of  the  notion  are  not  considered, 
much  more  information  is  available  than  if  only  such  parameters  as  the 
mean  or  the  variance  were  considered. 

Blackman  and  Tukey3  have  formulated  the  general  problem  of 
spectrum  analysis  of  digital  samples  and  provided  quantitative  guidance 
in  estimating  confidence  levels  of  the  analyses. 

In  recent  years,  cross-spectrum  analysis  has  teen  applied  exten¬ 
sively  in  several  branches  of  geophysics  and  has  proved  to  be  a  powerful 
tool  for  studying  certain  types  of  problems.  This  technique  has  teen 
effective  in  the  study  of  atmospheric  turbulence4  and  the  re  tractive- index 
structure  of  the  atmosphere.5  Munk,  Snodgrass,  and  Tucker6  have  applied 
cross-spectrum  techniques  extensively  in  angle-of'-arrival  and  coherence 
studies  of  ocean  waves. 

Description  of  the  methodology  of  spectrum  and  cross-spectrum 
analysis  and  its  application  is  widely  scattered  through  the  literature  of 
several  disciplines.  It  is  our  purpose  to  bring  together  in  one  document 


1  Superscript  miml>er>  identify  references  listed  ill  t>nd  of  report. 
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all  of  the  digital  formulation  required  by  most  engineers  and  applied 
physicists.  Elementary  derivations  of  the  basic  formulae  are  given.  The 
various  errors  that  may  contaminate  the  analysis  are  discussed  together 
with  the  confidence  limits  applicable  to  samples  of  finite  length.  Exam¬ 
ples  are  given  illustrating  the  application  and  the  step-by-step  procedure 
for  obtaining  spectra  and  cross  spectra.  Two  very  different  applications 
of  cross-spectrum  analysis  are  described,  and  the  significance  of  cross- 
power  and  coherence  is  discussed. 


POWER  SPECTRUM  ANALYSIS  •  GENERAL 


The  straightforward  way  to  obtain  the  power  for  a  continuous 
< noise)  process,  ,v( t ) ,  would  be  to  simply  Founer-analyze  the  data  sample. 
Thus 


E  ( co) 


1 

2  T* 


r  r  t* 

2 

‘  T* 

2 

/  ^x(t)  coscot  dt 

+ 

/  *xU)  sin  cot  dt 

l/-  T* 

T* 

— 

(1) 


where  the  total  sample  length  is  2  T*. 

If  we  digitize  v(/)  with  the  reading  interval  A  A  we  can  in  principle 
compute  7’*  A  /  harmonics  from  the  sample  and  the  bandwidth  is  therefore 
(2 T*)“‘.  As  T*  is  increased,  the  bandw  idth  correspondingly  decreases 
and  although  the  total  amount  of  data  is  increased,  the  information  content 
in  the  bandwidth  remains  the  same.  Therefore  the  scatter  of  spectral 
ordinates  is  not  reduced  and  we  find  the  surprising  result  that  increasing 
the  sample  length  does  not  improve  the  reliability  of  our  spectral  estimate. 
The  problem,  then,  is  to  find  a  way  of  getting  a  smooth,  reliable  estimate 
of  the  spectrum  of  the  population  and  to  be  able  to  state  confidence  limits 
for  the  estimates.  There  are  two  ways  of  doing  this.  We  might  carry  out 
the  harmonic  analysis,  computing  the  amplitude  of  the  individual  harmonies 
as  just  described,  and  then  smooth  the  spectrum  by  some  algebraic  smooth¬ 
ing  process.  This  would  be  the  equivalent  of  drawing  a  smooth  curve 
through  the  points.  Alternatively  we  might  carry  out  the  harmonic  analysis 
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of  the  autocorrelogram.  Tho  first  method  is  time-consuming  and  laborious, 
and  in  practice  the  best  method  is  usually  to  carry  out  a  harmonic  analysis 
of  the  autocorrelogram.  The  spectral  bandwidth  in  the  analysis  is  then 
(2-rm)  ~  where  t,„  is  the  maximum  time  lag  in  the  correlation  analysis. 

The  bandwidth  therefore  remains  constant  (for  a  given  total  lag)  and,  as 
the  sample  is  increased,  more  and  more  statistical  information  is  included 
in  each  frequency  band  leading  to  a  corresponding  improvement  in  the 
spectral  estimates. 


THE  COVARIANCE  FUNCTION 

The  power  spectrum  is  obtained  from  the  auto-covariance  function, 
and  the  cross-power  spectrum  from  the  cross-covariance  function.  For  two 
continuous  functions  x(t)  and  y(t)  whose  average  is  zero  the  covariance 
function  is  defined  to  be 

_  ,  rT* 

R  Y  (t)  -  x(t)  y(t  t  t)  T*=oa  /  x(t)  yW  +  t)  dt  (2) 

•VV  2  T*  J-T* 

where  2  T*  is  the  record  length,  the  bar  represents  the  mean,  and  t  is  the 
time  lag  of  x  relative  to  y.  The  function  R(t)  is  obtained  by  letting  t  vary 
from  0  to  sonic  maximum  value  tw  .  If  the  records  are  obtained  by  sampling 
at  an  interval  A  A  equation  2  is  approximated  by 


i  iY  —  p 


A  I 


x i i  t  p 


where  p  is  now  an  integral  number  of  A  t  steps  and  i  is  the  sequential  num¬ 
ber  in  the  sample  ranging  from  1  to  N. 

If  the  auto-covariance  is  computed,  xU)  and  y ( t)  are  identical  and 
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In  this  case  Rip)  is  symmetrical  around  p^O  and  only  one  side  of  the  cor¬ 
relation  function  need  be  calculated.  When  two  records  are  involved  both 
positive  and  negative  values  of  p  are  required,  and  p  must  take  all  integral 
values  0,  ±  1 ,  ±2 ,  ±m.  The  delay  time  t  is  given  by  p A  t  and  the  maxi¬ 

mum  delay  time  t,„  is  mA  t.  For  a  fixed  m,  jm  decreases  as  A  t  decreases. 


POWER-SPECTRUM  COMPUTATION 


For  a  continuous  covariance  function  Ri t)  the  power  spectrum 
E(co)  is  given  by  the  relation 


E  (<*>)=  —  f  Ri  t)  cos  COT  d  T  (5) 

t r  JO 

where  co  is  2  tt /'  and  f  is  cycles  per  unit  time. 

When  Ri t)  is  obtained  by  sampling  a  record  at  intervals  At,  the 
power  spectral  estimate  Eih)  is 


m 


I 

p  =  o 


Rip) 


(6) 


where 


8  = 


V2  when  p  -  0,  m 
1  when  0  <  p  <  m 


and  h  0,  1,  2,  ...  m.  The  frequency  in  cycles  per  unit  time  associated 
with  the  spectral  estimate  at  h  is  given  by 


f 


h 

2mA  t 


(7) 


The  frequency  bandwidth  A  f  in  cycles  for  each  h  is  f/h  or 
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(8) 


A  f  = 


1 

2mA  t 


Thus,  E if),  the  energy  in  a  bandwidth  of  one  cycle,  is  related  to  E(h)  as 

E(p  =  2mA  t  Eih)  (9) 


Systematic  Errors 


Two  systematic  errors  in  E ( //)  arise  because  of  truncation.  First, 
the  covariance  function  is  terminated  at  p  m.  Secondly,  the  sample  is 
truncated  because  it  is  finite. 

The  termination  of  Rip)  at  p  m  creates  “corner  effects”  in 
equation  5  and  produces  sidebands  in  E (/?).  This  is  partially  avoided  by 


forcing  R( t)  to  zero  at  p  -  m  by  multiplication  of  Rip)  in  equation  5  by 


(' 


*  cos 


The  other  truncation  error  arises  because  of  the  truncation  of  the 
sample  itself  and  is  a  result  of  the  finite  sample  length.  Thus  there  are 
artificial  abrupt  steps  in  practically  every  data  sample  subjected  to  cor¬ 
relation  analysis,  and  transients  arc  introduced  which  can  contaminate  the 
spectrum.  Prefiltering  of  the  records  to  remove  frequencies  lower  than 
those  being  analyzed  can  reduce  this  effect  but  it  cannot  be  entirely  elim¬ 
inated.  Furthermore  the  beginning  and  end  of  the  sample  should  be  chosen 
near  zero  values  of  the  variate  to  minimize  the  end  effects.  However 
truncation  is  removed  only  from  the  zero  order  lag  and  will  appear  in  sub¬ 
sequent  lags.  In  spectrum  analysis,  truncation  error  is  often  the  limiting 
factor  and  even  for  large  sample  sizes  it  will  usually  limit  the  validity  of 
the  spectrum  to  frequencies  whose  spectral  density  is  within  five  or  six 
orders  of  magnitude  of  the  maximum  spectral  density  in  the  spectrum. 

If  the  spectral  density  varies  rapidly  with  frequency,  the  energy 
contributed  to  the  spectrum  by  the  side-lobes  of  the  spectral  window  on 
the  low-energy  side  of  the  center  frequency  will  lx*  greater  relative 
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-Vf  =  *i  +  b  i 


where  b  lies  between  +1  and  -1.  Consider  one  Fourier  component. 
X;  =  C  cos  i.  Then 


.X  j  _  j  =  C  cos  ^  TAJ  ( / _  1 ) 
and  y,  represents  a  sum  of  the  form 


(10) 


Ill) 


-i+aj  ^  b  cos  f—^ i -  p^ 


(12) 


where  a  -  1 ,  a  --  0,  and  p  =  - 


2  ttA  t 
T 


.  The  amplitude  of  y  is  the  resultant 


of  the  sum  of  two  vectors  of  magnitude  a ,  b  and  phases  a,  (3  so 

j  -  a1  t  6J  +  2a6  cos  (a  -(3 ) 

Therefore  the  filter  response  of  this  prewhitening  scheme  is 


(13) 


F 2  -  1  -+  b2  ■  2  b  cos  &  t 

T 


(14) 


to  the  true  spectrum  than  that  contributed  by  the  lobe  on  the  high-energy 
side.  The  effect  is  to  shift  energy  in  the  true  spectrum  toward  lower-energy 
regions  in  the  spectrum.  It  may  therefore  be  advisable  to  arrange  the  sam¬ 
ple  so  that  the  energy  is  more  or  less  evenly  distributed  through  the 
spectrum  before  carrying  out  the  actual  analysis.  This  is  sometimes 
called  “prewhitening.”  One  method  of  prewhitening6  modifies  the  data 
sample  v.  in  the  following  way.  Let 


If  T  x,  F  =  1  i  5.  If  T  :  2  A  /  (the  highest  frequency  in  the  analysis), 
FI-  b.  Therefore,  if  b  -  -0.75,  the  energy  at  high  frequencies  is  in¬ 
creased  by  a  factor  of  40  lb  and  the  energy  at  low  frequencies  is  decreased 
to  1  10  its  actual  value. 
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The  effects  of  averaging  within  each  data  point  must  also  be  con¬ 
sidered.  The  time  series  is  devoid  of  the  highest  frequencies  if  the  obser¬ 
vations  are  from  averaged  data.  In  other  words,  averaging  acts  like  a  low- 
pass  filter.  The  nature  of  this  filter  function  is  well  known.  If  E(/>  is  the 
spectrum  of  some  variable  .v,  the  spectrum  of  the  mean  .v  averaged  over  a 
period  ta  is  given  by 


Ky</-> 


(15) 


Therefore  it  is  easy  to  obtain  the  spectrum  of  x  if  the  spectrum  of  .v  is 
known. 

Tukey  coined  the  phrase  “aliasing"  for  the  phenomenon  that 
causes  energy  from  one  frequency  to  appear  under  the  "alias"  of  another 
frequency.  The  effect  is  due  to  the  discrete  sampling  or  digitizing  of 
records  for  analysis  by  digital  computer.  The  highest  frequency  in  the 
analysis  has  a  period  twice  that  of  the  sampling  interval.  However,  if 
higher  frequencies  are  present  in  the  record,  these  higher  frequencies  may 
still  appear  in  the  analysis  but  disguised  as  a  lower  frequency.  Consider 
the  following  figure: 


wwww 


Suppose  that  this  sinusoidal  record  of  frequency  f  is  sampled  at  the  points 
circled.  Then  the  highest  frequency  in  the  spectrum  analysis,  (2  A  0"  l, 
will  be  about  half  the  frequency  prominent  in  the  record.  However,  since 
tlx.'  sampling  interval  A  /  is  a  little  less  than  the  wavelength  in  the  record, 
the  sampling  points  fall  first  on  the  crests,  later  in  the  troughs,  then  again 
on  the  crests  pnxlucing  the  dashed  curve  of  frequency  f  -  ( A  t)~l.  The 
original  record  contained  no  such  frequency.  Instead,  it  was  intnxluced 
artificially  by  the  large  spacing  between  observations.  Frequencies  higher 
than  those  analyzed  appear  under  the  alias  of  the  lower  frequency 
/  (A  0~\  and  resemble  a  modulation  not  unlike  a  “beating"  of  the  fre¬ 
quencies  /  and  (A  /)"*.  The  effect  of  aliasing  is  most  pronounced  at  tin* 


high-frequency  end  of  the  spectrum.  Therefore,  one  method  of  minimizing 
the  effects  of  aliasing  is  to  “over-sample”  the  record  and  then  ignore  the 
higher  spectral  ordinates. 


Random  Errors 

Random  errors  due  to  finite  sampling  must  also  be  considered. 
Different  samples  from  the  population  would  lead  to  different  spectral 
estimates,  each  representing  a  summation  over  the  harmonics  in  the  sam¬ 
ple.  If  samples  are  drawn  at  random  from  a  statistically  stationary  record 
and  if  the  variate  is  normally  distributed,  the  spectral  estimate  for  a  given 
frequency  will  vary  about  the  population  spectrum  according  to  the  chi- 
square  distribution  which  is  widely  tabulated.  The  value  of  chi-square  for 
a  given  significance  depends  only  on  the  number  of  degrees  of  freedom, 
which  is  defined  as  the  number  of  independent  normally  distributed  quan¬ 
tities  that  are  summed  -  -in  our  case  the  number  of  harmonics  in  the  sum¬ 
mation.  The  effective  width  of  the  main  lobe  of  the  window  of  the  Tukey 
filter  (Hanning  function)  is  2.  Since  this  will  pass  ,V  m  harmonics  through 
each  of  the  two  windows,  the  number  of  degrees  of  freedom  is  approximately 
2 i\v'm.  Actually,  Tukey  suggests  that  the  number  of  degrees  of  freedom, 
v,  is  more  like  2.5(,\Vm)-l/2  but  he  recommends  the  use  of  2(\/m) -1/2 
because  of  uncertainties  such  as  the  assumption  of  a  normal  distribution 
of  the  variate,  the  appropriateness  of  the  chi-square  distribution,  and  the 
statistical  stationarity  of  the  record.  The  confidence  limits  for  the  spec¬ 
tral  ordinates  are  related  to  the  degrees  of  freedom  in  figure  1. 

It  is  evident  that  the  confidence  limits  depend  directly  on  sample 
size,  /V,  and  on  the  total  amount  the  record  is  lagged,  m.  The  expression 
2  <  S/m)  - 1/2  relating  degrees  of  freedom  to  sample  size  and  total  lag  is 
important,  and  several  things  can  bo  immediately  pointed  out.  First  of  all. 
when  the  sample  size  is  large  compared  to  the  total  lag,  v  is  nearly  pro¬ 
portional  to  the  ratio  of  sample  size.  ,\,  to  maximum  lag  length,  m.  There¬ 
fore  this  ratio  is  the  important  factor  in  determining  the  confidence  limits 
of  the  spectrum  analysis.  Furthermore,  the  maximum  lag,  t,„.  determines 
the  lowest  frequency  (longest  period)  that  can  be  analyzed  in  the  data 
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Figure  1.  Ninety  percent  confidence  limit*  relnt ing  degrees  of  freedom 
(u  the  ratio  of  computed  energy  to  irne  ciiergv  density.  Nine  times  out  of 
ten  the  computed  >pc<  tl'al  delis  1 1  \  wilt  lie  between  the  curves  when  ex¬ 
pressed  us  a  ratio  to  the  correct  population  density.  tAftci  Blackman 
and  Tuke\  ;  see  ref.  b  * 
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sample  since  /’mjn  -  (2tot)  —1,  while  the  digitizing  interval,  At,  determines 
the  highest  frequency  since  fmax  =  (2 A  f)~‘-  This  means  that  the  lowest 
frequency  of  interest  in  the  analysis  determines  the  sample  length,  2 T*, 
required  in  order  to  obtain  a  satisfactory  confidence  level  in  the  spectrum 
analysis.  This  is  a  point  to  be  emphasized  because  sometimes  attempts 
are  made  to  get  increased  reliability  by  just  sampling  the  record  more  often 
to  get  a  larger  \  without  increasing  T*.  The  sampling  frequency  determines 
the  highest  frequency  capable  of  detection  in  the  record  just  as  the  maxi¬ 
mum  lag  (time)  determines  the  lowest  frequency.  Since  the  highest  fre¬ 
quency  in  the  spectrum  will  be  (2  A  d-1,  increasing  the  sampling  rate 
(that  is  decreasing  the  sampling  interval  A  t)  simply  extends  the  analysis 
to  higher  frequencies  and  the  confidence  limits  of  the  spectrum  remain  un¬ 
changed  if  T*  Tm  is  unchanged. 


Stationarity  and  Filtering 

The  foregoing  discussion  has  pointed  out  the  adverse  effects  pro¬ 
duced  by  frequencies  outside  the  range  of  those  being  considered  in  the 
spectrum  analysis.  They  lead  to  errors  from  aliasing  and  to  erroneous 
confidence  estimates  related  to  nonstationarity. 

A  sample  must  be  stationary  if  it  is  to  be  considered  representative 
of  the  population.  Stationarity  normally  means  that  the  lower-order  statis¬ 
tics  should  not  vary  throughout  the  record.  If  the  mean  and  variance  of 
the  record  are  constant  everywhere  in  the  record  it  is  generally  considered 
stationary.  Consider  the  record  fE  shown  at  the  top  of  figure  2.  fE  is  an 
average  of  fE  about  local  noon  and  indicates  the  average  ionospheric 
electron  density.*  The  record  gives  fE  for  180  days  and  contains  a  grad¬ 
ual  decrease  in  the  running  average  (averaged  over  five  to  ten  days)  which 
is  caused  by  the  annual  change  in  the  solar  zenith  angle.  Thus  the  record 
of  fE  is  not  stationary  relative  to  the  mean,  the  first-order  statistic. 


*  fE  is  the  critical  frequency  for  radio  penetrntion  of  the  R  Region  of  the 
ion  os  phe  re. 
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Figure  2.  A  portion  of  a  record  of  fE  (see  text  for  definition!.  The 
upper  record  is  the  unfiltcred  raw  data:  the  lower  record  is  data  filtered 
using  equation  A-7  with  m  --  100  days. 
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If  the  record  of  fE  extended  over  four  years,  four  annual  cycles 
would  be  present.  The  record  should  not  be  considered  stationary  because 
!0  to  20  cycles  at  the  lowest  frequency  in  a  record  is  considered  neces¬ 
sary  for  stationaritv. 

An  eleven-year  record  of  fE  would  show  a  longer-period  variation 
associated  with  the  1  l-to-12-year  sunspot  cycle,  and  again  the  record 
would  l)e  nonstationary.  Further  extension  of  the  record  might  reveal 
longer-period  variations  which  might  be  associated  with  climatic  changes. 
Most  geophysical  data  are  correspondingly  nonstationary  and  should  be 
converted  to  a  stationary  form  relative  to  the  mean.  Nonstationarity  in 
the  mean  can  be  removed  by  numerically  prefiltcring  the  data  as  shown  by 
the  bottom  record  of  figure  2. 

Sonic  data  show  a  change  of  variance  throughout  the  record  and 
the  record  is  considered  nonstationary.  Nonstationarity  in  the  variance 
cannot  be  removed  by  prefiltering.  Unless  the  variation  is  judged  to  be 
excessive  by  some  method,  the  nonstationarity  of  the  variance  is  usually 
ignored.  Nonstationarity  in  the  variance  or  other  higher-order  statistic- 
will  not  Ik?  considered  here. 

Low-  or  high-frequency  components  in  a  record  may  be*  removed  by 
mechanical-electrical  devices  during  recording  or  by  mathematical  methods 
after  recording.  Inductors  or  capacitors  arc  often  used  to  remove'  low  fre¬ 
quencies  and  a  time  constant  may  be  introduced  by  an  RU  circuit,  to  remove* 
high  frequencies  prior  to  recording  a  signal.  Various  mathematical  filters 
are  available  for  operating  on  recorded  data,  depending  on  what  frequencies 
one*  desires  tee  remove  or  reduce  in  amplitude.  The  important  requirement 
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is  that  the  frequency  response  of  the  filter  must  lx?  known  so  that  the 
analysis  of  the  spectral  information  can  be  objective.  Mechanical  or  elec¬ 
trical  filter  responses  will  not  be  examined.  Examples  of  two  convenient 
numerical  filters  are  given  in  Appendix  A. 

The  foregoing  discussion  of  the  main  elements  involved  in  auto¬ 
correlation  and  power-spectrum  analysis  leads  now  to  consideration  of  a 
simple  example. 


Sample  Problem  (Power  Spectrum] 


Consider  the  record  of  geomagnetic  fluctuations  shown  in  figure 
3.  The  first  step  in  any  spectrum  analysis  is  to  decide  the  range  of  fre¬ 
quencies  of  interest.  In  the  present  record,  it  is  obvious  that  we  have  a 
strong  60-Hz  component.  We  will  therefore  wish  to  terminate  the  spectrum 
analysis  at  a  frequency  much  less  than  60  Hz  --  say  30  Hz.  In  this  data 
sample,  we  are  looking  for  earth-cavity  resonance  effects.  Theory  offers 
some  guidance  and  we  know  that  the  fundamental  frequency  should  be 
some  few  Hz.  We  might  therefore  decide  to  analyze  the  spectrum  down  to 
a  frequency  of  perhaps  0.5  Hz.  This  frequency  determines  the  value  of 
to  be  used  in  the  analysis.  t,„  not  only  determines  the  longest  period 
that  we  can  analyze  in  the  record,  but  it  also  determines  the  reso¬ 
lution  in  the  spectrum  analysis  (A  f  =  (2-rm)~l).  The  smaller  the  maximum 
lag,  m,  compared  to  the  number  of  observations,  .V,  the  more  smoothing  is 
obtained  in  the  spectrum.  Also,  it  is  important  to  remember  that  less  com¬ 
puter  time  is  required.  However,  a  small  value  of  m  limits  the  resolution 
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Figure  3.  Portion  of  magnetometer  record  showing  micropulsations  of 
earth  s  magnetic  field.  Record  is  contaminated  with  f50-cyclo  ''noise'  . 
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and  may  obscure  significant  characteristics  in  the  spectrum  of  the  popu¬ 
lation.  If  we  wish  to  examine  frequencies  as  low  as  0.5  Hz,  the  maximum 
lag,  m,  in  the  correlation  analysis  must  extend  over  at  least  a  1 -second 
time  interval.  At  the  other  end  of  the  spectrum,  if  we  wish  to  examine 
frequencies  as  high  as  30  Hz,  we  must  sample  the  records  at  twice  that 
frequency,  or  60  samples  per  second.  At  a  sampling  rate  of  60  Hz,  there 
are  60  lags  in  a  time  interval  of  1  second,  therefore  m  60.  We  have  now 
selected  the  range  of  frequencies  of  interest  and  the  maximum  lag  number 
m  in  the  correlation  analysis.  It  remains  to  determine  how  long  a  data 
sample  will  be  required  to  achieve  the  ordinate  confidence  we  require. 
Suppose  we  decide  that  100  degrees  of  freedom  is  adequate.  This  means 
that  90  percent  of  the  time  the  estimated  ordinate  value  will  lie  between 
approximately  0.77  and  1.25  of  the  true  spectral  value  for  the  population. 
For  100  degrees  of  freedom  and  a  lag  of  60,  we  find  that  our  required  sam¬ 
ple  length  is  3000.  We  are  now  ready  to  carry  out  the  spectrum  analysis. 
We  will  usually  wish  to  express  the  spectral  values  as  the  variance  of  the 
data  per  cycle  per  unit  time.  That  is,  if  the  record  is  recorded  in  seconds 
we  will  usually  wish  to  express  the  energy  as  variance  per  cycle  per  sec¬ 
ond.  This  has  been  done  in  the  sample  shown  and  the  resulting  units 
are  < m i  1 1  i  y>J/Hz. 

In  geophysical  problems,  it  often  happens  that  the  spectrum  falls 
off  rapidly  as  we  proceed  toward  higher  frequencies,  and  it  may  in  fact 
drop  several  orders  of  magnitude  within  the  range  of  frequencies  analyzed. 
It  is  therefore  often  convenient  to  plot  log  E  as  function  of  frequency. 
Alternatively,  it  may  be  convenient  to  multiply  each  spectral  value  by  its 
frequency  to  obtain  a  spectral  function  with  dimensions  of  variance  and  a 
reduced  range  of  variation.  The  spectrum  can  usually  then  be  conveniently 
plotted  on  linear  paper  without  the  undue  suppression  of  spectral  lines  that 
may  occur  at  the  higher  frequencies.  This  is  the  method  chosen  for  pres¬ 
entation  in  figure  4. 
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CROSS  SPECTRUM  COMPUTATION 


The  power  spectrum  provides  information  about  the  frequency 
content  of  a  time  series.  Cross-spectrum  analysis  provides  information 
about  the  relationship  between  two  or  more  time  series  in  the  frequency 
domain,  i.e.,  the  coherence  and  phase  lag.  Goodman7'8  discusses  the 
theory  and  principles  of  cross-spectrum  analysis.  The  multiple  coherence 
is  analogous  to  the  multiple-correlation  coefficient.  In  a  two-parameter 
analysis  the  phase  lag  is  analogous  to  the  intercepts  of  the  regression 
line  of  best  fit  in  the  correlation  analysis.  Just  as  significance  levels 
are  required  in  a  correlation  analysis,  confidence  levels  are  also  required 
in  cross-spectrum  analysis.  This  section  will  present  some  elementary 
concepts  of  coherence  and  phase  lag  and  present  some  examples  in  which 
confidence  levels  are  used. 

Just  as  the  power  spectrum  can  be  obtained  by  Fourier-transform- 
ing  the  auto-covariance  function,  the  cross  spectrum  is  obtained  by  trans¬ 
forming  the  cross-covariance  function.  If  the  cross-covariance  function 
is  asymmetrical,  both  even  and  odd  terms  are  required  to  represent  it  and 
the  cross  spectrum  will  therefore  be  complex.  Therefore 
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Now  /?,,( t)  can  lx*  decomposed  into  a  symmetric  term  and  an  antisymmetric 
term  as 
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where  the  symmetric  and  antisymmetric  terms  are  as  shown  schematically: 
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The  smoothed  spectral  estimates  of  C(/>  and  Q(f)  at  /i  (corres¬ 
ponding  to  frequency  /  =  h/(2mAt»  can  be  numerically  calculated  from 
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When  records  1  and  2  are  identical  Ru  -  RM  =  R22  is  symmetrical 
about  t-  0,  and  C<P  becomes  the  power-spectral  estimate  while 
vanishes. 

The  phase  relationship  between  frequencies  in  records  being  cross¬ 
spectrum  analyzed  is  given  by 


while  the  cross-power  is 
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The  coherence  is  defined  as  the  square  of  the  cross-power  divided 
by  the  product  of  the  powers  in  the  two  records  being  analyzed,  or 


Coherence  measures  the  fraction  of  the  variance  of  one  variable 
(record  1)  that  can  be  specified  by  another  variable  (record  2)  because  of 
phase  coherence  due  to  some  direct  or  indirect  association  between  the 
records  at  frequency  co.  The  coherence  should  vary  only  from  0  to  1,  but 
coherence  quickly  becomes  unreliable  if  the  degrees  of  freedom  in  the 
analysis  arc  insufficient  and  may  oven  exceed  unity  if  the  confidence  in 
the  numerical  analysis  is  low. 

When  more  than  two  recot ds  are  involved,  multiple  coherence  can 
be  calculated  with  interpretations  similar  to  multiple-correlation  coeffi¬ 
cients.  Multiple  coherence  will  not  be  considered  here. 

There  is  a  subtle  difficulty  that  appears  when  wo  work  with  data 
of  finite  length.  The  problem  is  very  well  stated  by  Madden,"  as  follows: 

“The  cross-correlation  of  two  transients  does  not 
lose  any  of  the  power  of  the  transients,  and  this  would  ap¬ 
pear  to  guarantee  a  coherency  of  1 .  This  difficulty  is  over¬ 
come  by  making  the  spectral  window  average  together  the 
cross-power  of  a  group  of  neighboring  frequencies.  When 
the  two  noise  signals  are  random  to  each  other,  the  phase's 
of  each  frequency  component  are  random,  and  a  vector 
addition  of  the  cross-powers  of  a  group  of  neighboring  fre¬ 
quencies  tends  to  cancel.  When  the  two  noise  signals  arc 
coherent,  the  phases  of  each  frequency  in  the  cross-power 
are  those  of  the  phase  spectrum  of  the  linear  operator 
which  generates  one  signal  from  the  other.  If  this  spec¬ 
trum  is  slowly  varying  as  a  junction  of  frequency,  the  ad¬ 
dition  of  neighboring  cross-powers  docs  not  cancel  out. 

This  necessitates  some  rare  in  evaluating  coherency  esti¬ 
mates  if  time  delays  are  involved,  or  if  very  sharp  frequency 
structure  exists  in  the  relationship  between  the  data.” 

The  probability  of' obtaining  a  particular  coherence  by  chance 
(the  level  of  confidence*  must  be  sjx’cified  if  the  statistical  analysis  is 


to  be  considered  complete.  The  confidence  level  implies  a  reliability 
level  for  coherence.  The  probability  density  P  of  coherence  has  been 
considered  by  Goodman7  and  is  approximately  given  by: 
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where  coh*  is  the  population  coherence.  Alexander  and  Vok,c  have  tabu¬ 
lated  the  cumulative  distribution  of  sample  coherence  as  a  function  of  the 
degrees  of  freedom  [y  =  (2.V  m)  -0.5]  and  the  true  or  population  coherence. 
Unfortunately,  v  ranges  from  2  to  only  21  in  the  tables.  Population  coher¬ 
ence  ranges  from  0  to  0.99.  The  true  coherence  is  the  coherence  expected 
from  the  population  having  v  degrees  of  freedom.  Figure  5  shows  the  cumu¬ 
lative  distribution  function  of  coherence  for  5,  10,  15,  and  20  degrees  of 
freedom  when  the  population  coherence  is  zero.  Figure  6  shows  the  cumu¬ 
lative  distribution  function  of  coherence  for  population  coherences  of  0,  5. 

10,  15,  20,  30,  and  40  when  v  is  10.  These  data  were  taken  from  reference  10. 

The  probability  P  that  phase  will  lie  within  a  certain  increment 
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Hu1  limits  for  P  -  0.95  are  shown  in  figure  7  (taken  from  ref.  b> 


Figure  7.  Nincly-fi vi*  percent  confidence  limits  for  phase  angles 
as  functions  of  degrees  of  freedom,  v  .  Thus  for  v  =  100  degrees 
ol  trec'dom  there  are  nineteen  chances  in  twenty  if  coh*  =  0.5  that 
the  true  phase  is  within  0  25  radian  of  ihe  computed  phase. 

(From  Munk  cl  a!.:  see  ref.  (>.' 


Sample  Problem  Applying  Cross-Spectrum  Analysis  to 
Time-Related  Functions 


The  electron  density  in  the  ionosphere  is  known  to  be  controlled 
primarily  by  the  intensity  of  the  solar  ionizing  energy.  Aside  from  the 
diurnal  and  annual  variation,  the  electron  density  varies  with  the  sunspot 
number  which  has  a  long-term  II -year  cycle  and  a  short-term  27-day  cycle 
A  1 5-day  or  semilunar  cycle  has  also  been  found.  A  study  by  Noonkester 
was  made  to  determine  the  behavior  of  two  ionospheric  layers  < F  and  F> 
at  the  15-  and  27-day  |x>ri<xls.  Daily  measures  of  [F  and  fF  represent  the 
variation  of  electron  density  at  the  F  and  F  levels.  They  were  obtained 
for  2722  consecutive  days  and  the  annual  cycle  was  removed  (see  filtered 
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data  in  the  sample  record,  fig.  2).  The  power  spectrum  and  the  coherence 
spectrum  were  then  obtained  for  fE  and  fF. 

The  power  spectra  of  j  E  and  fF  showed  a  significant  maximum 
near  a  period  of  27  days,  but  did  not  show  any  maximum  near  the  15-day 
period.  The  power  spectra  of  each  quarter  of  the  data  showed  that  the  27- 
day  period  was  dependent  on  the  portion  of  the  records  considered. 

A  coherence  spectrum  between  fE  and  fF  (fig.  8)  showed  that 
they  were  related  near  the  27-  and  15-day  periods.  The  spectral  estimates 
have  42  degrees  of  freedom  and  are  separated  by  4x10“"'  cycles  day.  The 
confidence  limits  cannot  be  determined  from  the  tables  of  Alexander  and 
Yok10  because  these  tables  do  not  consider  v  greater  than  21.  However,  an 
estimate  of  the  significance  level  was  made  at  a  v  of  21.  The  population 
coherence  was  assumed  to  be  the  average  of  the  125  estimates  of  coher¬ 
ence  in  the  spectrum.  According  to  the  tables,  at  21  degrees  of  freedom 
a  coherence  of  0.35  or  greater  has  a  probability  of  less  than  0.01  of  occur¬ 
ring  by  chance.  Certainly  the  coherence  estimates  in  figure  8  at  the  27- 
and  15-day  periods  have  a  small  probability  of  occurring  by  chance. 
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The  value  of  the  coherence  is  shown  in  the  above  example.  The 
power  spectra  of  fE  and  fF  indicated  that  the  amplitude  of  the  15-day 
period  was  at  the  background  noise  level,  but  the  coherence  between  these 
variables  indicated  that  they  were  related  at  the  15-day  period.  Although 
the  amplitude  was  small  at  the  15-day  period,  fE  and  fF  varied  coherently 
at  the  15-day  period.  This  is  equivalent  to  using  phase  coherence  to  ob¬ 
tain  a  signal  out  of  a  noisy  record. 

According  to  equations  22  and  26,  fE  was  found  to  precede  fF  by 
V.2  days  near  the  27-day  period  within  confidence  limits  of  ±  2.0  days  90 
percent  of  the  time.  fE  and  fF  were  found  to  be  in  phase  at  the  15-day 
period  and  within  confidence  limits  of  ±1.2  days  90  percent  of  the  time. 


Application  cf  Cross-Spectrum  Analysis  to 
Spatially  Separated  Recordings 

In  this  application,  we  consider  two  separate  records  of  a  time- 
varying  function  received  at  two  spatially  separated  stations  A,,  A2  shown 
schematically  in  figure  9.  The  function  might  be  generated  by  the  steady 
drift  across  the  stations  of  a  spatially  varying  medium  perturbed  by  turbu¬ 
lent  processes;  or  it  might  be  generated  by  a  wave  motion  of  significant 
bandwidth  propagating  through  the  medium  where  the  different  frequency 
components  may  propagate  in  different  directions  with  different  velocities. 
Figure  9A  shows  the  case  for  a  single  Fourier  component  of  a  wave  struc¬ 
ture  propagating  in  the  x  direction.  In  figure  9A  let  t0  be  the  time  lag 
between  the  arrival  of  the  wave  crest  at  A,  and  A2.  From  equations  22 
and  23 


C(  co)  =  |i 

fl/j  2  ^  CoH  COS  CO  To 

(27) 

Q(co>  =  1 

Eis  (co)|sin  co  To 

(28) 

If  T0  =  0,  Q(f)  -  0.  If  t„  is  not  zero,  tan  cot,,  -  Q( co)/C(co)  gives  the 
phase  lag  between  the  records  corresponding  to  the  distance  of  travel 
A’t0  L  sin  0.  Since  co  (2tt'A)V  kV,  tan  i  (2tt.  "K)(L  sin  0)1  Q(0  (’(/"), 
so  the  wavelength,  A,  and  the  average  angle-of-nrrival,  0  (see  fig.  9A>,  can 
be  found  as  a  function  of  frequency  f. 


A.  GEOMETRY  OF  WAVE  TRAIN  OF  WAVELENGTH.  A.  PASSING 
RECORDERS  LOCATED  AT  A,  AND  A2  SEPARATED  BY  A 
DISTANCE  L.  THE  LINE  JOINING  THE  STATIONS  MAKES  THE 
ANGLE  0  WITH  THE  v  AXIS.  9  IS  THEREFORE  THE  AVERAGE 
ANGLE  OF  PROPAGATION  OF  THE  WAVE  TRAIN  WITH  RESPECT 
TO  THE  ORIENTATION  OF  THE  RECORDER  PAIRS. 
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B.  GEOMETRY  OF  TWO  INTERFERING  WAVE  TRAINS  WHOSE 
NORMALS  MAKE  THE  ANGLE  9 WITH  THE  a  AXIS.  THEY  NUKE 
ANGLES  6  ±  ©  WITH  RESPECT  TO  THE  NORMAL  TO  THE  LINE 
OF  STATION  SEPARATION. 
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The  cospectrum  C(f)  corresponds  to  the  in-phase  frequencies  of 
the  spectrum,  and  the  quadrature  spectrum,  Q(f),  corresponds  to  the  out- 
of-phase  components.  Since  the  cosine  transform  (cospectrum)  is  an  even 
function,  the  corresponding  correlation  function  must  he  symmetric.  The 
correlation  function  corresponding  to  the  sine  transform  (quadrature  spec¬ 
trum'  is  anti-symmetric.  Since  the  autocorrelation  function  is  always 
symmetric  for  a  stationary  time  series,  the  corresponding  quadrature  spec¬ 
trum  is  always  zero,  and  the  cospectrum  is  identical  to  the  power  spectrum. 
A  cross-correlation  function  is,  in  general,  nonsymmetric  about  0.  and 
Q(/‘)  is  not  zero.  The  argument  cot0  is  then  the  phase  lag  between  the  two 
records  which  were  cross-correlated  and  represents  the  phase  cross-spec¬ 
trum  between  the  records. 

The  normalized  magnitude  squared  of  the  cross-spectrum  between 
two  records  is  called  the  coherence.  It  is  sometimes  stated  that  coherence 
is  like  a  correlation  function  representative  of  a  frequency  band  in  the 
spectrum.  This  can  be  a  useful  concept  if  used  cautiously,  but  factors 
that  can  reduce  the  coherence  between  records  are  now  considered. 

Suppose  we  are  receiving  some  signal  at  two  recording  sites 
spatially  separated.  Suppose  the  signal  is  modulated  by  temporal  changes 
in  the  propagation  medium  and  we  wish  to  obtain  information  about  these 
changes  by  examining  the  coherence  of  the  fluctuations  in  the  two  records. 
The  coherence  will  be  degraded  by: 

1.  Random  or  noisciike  phenomena  such  as  turbulence  in  the 
neighborhood  of  each  recorder. 

2.  Anything  that  causes  the  same  frequency  to  have  a  variable 
travel  time  between  the  recorders  such  as  a  changing  direc¬ 
tion  of  propagation  during  the  sampling,  or  nonlinear  processes 
in  wave  phenomena. 

3.  Source  or  sink  located  between  recorders. 

-1.  The  same  frequencies  coming  simultaneously  from  different 
directions  as  would  lie  the  case  if  wavelike  phenomena  in 
the  medium  had  a  significant  "lieam  width  "  or  a  distributed 
source.  Then  interfere  nee  Ix’twocn  waves  coming  from 


different  directions  at  the  same  frequency  will  “smear”  the 
phase  difference  between  recording  stations  and  reduce  the 
coherence. 

The  last  mechanism  is  of  both  theoretical  and  practical  interest  and  has 
been  used  by  oceanographers  (i.e.,  Cox;  see  ref.  12)  to  study  the  “beam 
width”  of  ocean  waves.  The  problem  will  now  be  formulated. 

It  is  clear  from  the  preceding  development  of  the  expressions  for 
the  cross-spectra  that  the  \J coherence  is  just  the  normalized  modulus  of 
the  complex  cross-spectrum  and  the  phase  angle,  cot0,  is  its  argument. 

Suppose  now  that  waves  of  the  same  frequency  arrive  at  the  re¬ 
corders  from  different  directions,  0  +  9,  where  0  is  the  average  direction 
(see  fig.  9B).  Let  the  co-  and  quadrature  spectra  represent  the  average 
of  spectral  components  over  all  directions  of  arrival.  Then  cot0  =  kL  sin 
(0  +  9)  =  /;L(sin  0  cos  9+  cos  0  sin  9).  Let  x  and  y  be  the  components  of 
L  in,  and  normal  to,  the  average  direction  of  propagation,  respectively. 
Then  from  equations  27  and  28  we  have 


C(„,  =  _»_/** 

2A9./-A9 

1  rA  $ 
Q(co)  =  / 

2A9J-A9 


E,  j  ( co,  9)  cos  ( kx  cos  9  +  ky  sin  9)  t/  9 


Eu  (00,9)  sin  ( kx  cos  9+  ky  sin  9)0/9 


(29) 

(30) 


Now  if  we  consider  the  coherence  to  be  degraded  only  by  the  interference 
of  waves  of  the  same  frequency  coming  from  different  directions  with  dif¬ 
fering  co  To,  equations  27  and  28  show  thntlEj2  =  E,  Ea.  Furthermore,  if 
the  wave  structure  is  spatially  homogeneous,  E,  (co,  9)  =  Ea  (co,  9)  = 

E(co,  9).  If  E  (00,9)  is  assumed  to  be  uniform  for  -A  9^  9-  +  A 9 and  zero 
elsewhere  it  can  be  removed  from  the  integral  so  that  equations  29  and  30 
give: 

(C  ♦  iQ)  (co,  A 9)  _  _J —  f  exp  l-//.’(.v  COS9+ v  sin9>U/9  (31) 
E  (co)  2A9J--A9 

'finis  coherence,  equation  21,  which  is  just  the  square  of  the  cross-power 
normalized  to  the  auto  powers,  can  lx'  used  as  a  measure  of  the  “beam 
width"  of  wave  components. 
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An  alternative  concept  of  the  dependence  of  the  cross-power  on 
wave  pattern  is  to  define  a  “radiation  pattern”  F(q>)  of  the  recorder  array. 
Stated  formally 


so  that 


F(q>)  J E(co,q>)  d<p  =  J E(  co,  9)  exp  [~ik(: <  cos  9  +  ysin<p)]d(p 
Je( co,  9)  exp  [~-ik(x  cos  9+  y  sin 9)]  dy 


F(  9)  = 


j^Ei go,  9)  d<p 


(32) 


(33) 


where  E  (00,9)  is  the  directional  spectrum  of  the  radiation. 

If,  as  before,  E< co)  is  constant  for  -A9<  9<  A 9 and  zero  else¬ 
where,  E(oo)  can  be  removed  from  the  integral  and  the  limits  will  be  -  A  9 
and  A  9  so  that 


1  f 


1  rA  9 

F(A9>=  — —  J  ^  exp  [~ik(x  COS9+  y  sin <p)]  dtp 


(34) 


and  we  thus  identify  FU\<f)\  with  \f coherence  by  equation  31.  Equation 
34  can  be  written 

1  rAcp 

F(A9)  -  - -  /  [cos (k\  sin9>  —  /  sin  ( A? v  sin®)] exp(-/7j.vcos 9) da> 

2  A  9  ■/-  A  9  •  '  (35) 


But  9  was  chosen  as  the  angle  off  the  x  axis  so  the  integral  of  the  sin 
( hy  sin  9)  function  from  -A 9 to  A 9  is  zero. 

Therefore 

1  rA<F 

F(A9)  =  — /  cos  (fry sin9)  cxp(  —  ikx  cos 9)  d 9  (36) 

A9^0 

A  limiting  case  that  is  interesting  to  consider  is  that  of  a  very 
narrow  beam  (A 9--.  ■  1.0).  Then  9 is  small  and  equation  36  is  approximately 
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=  2i l!y<fexp  [-'**  ('  -ty 


r/< p 


=  £^r. /Xb  -Vf>]  ^ 


exp(— ikx) 
4  A  9 


rA»  vn|/tvr/  vy  ivv 

-A»exp  T[r‘)" (“) 


d  9 


(37) 


r/  9 


Transforming  variables  so  that  r+3  (2x  AK9  +  V  x)3  and  being  careful  in 

rearranging  terms  according  to  appropriate  limits,  equation  37  can  be 
written  as  follows: 


Ft  A  9) 


i'2dv 


j: 


cos  -  r2  dr 
9 


(38) 


3  dr  +  i  /  sin  ^  r2  dr 
Jo  2 


which  is  the  convenient  form  since  the  integrals  are  the  Fresnel  integral 
which  is  commonly  tabulated.  To  arrive  at  this  form,  note  that  both  the 
cosine  and  sine  terms  are  even  functions. 

For  station  separation  in  the  x  direction  only  (stations  located 
along  the  direction  of  propagation)  y  is  zero  and  the  equation  simplifies 
to  just  two  Fresnel  integrals.  The  result  is  plotted  as  the  black  curves 
corresponding  to  various  beam  widths,  A 9,  in  figure  10.  Note  that  the 
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left-hand  scale  is  to  be  used  and  its  great  expansion  shows  immediately 
that  coherence  falls  off  very  slowly  for  separations  along  the  direction  of 
propagation  even  for  fairly  large  beam  widths. 


1.00 
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!'  inure  10.  Dependence  <>f\  coherence  <>n  recorder  sc j>; i r: »( ion  and 
wavelength,  v  is  separation  in  the  direetion  of  propagation.  \  is 
ion  porpoixl  iruliir  to  tin*  dimtion  nf  proposition. 


For  station  separation  in  the  y  direction,  \  0.  (The  line*  of 

st.it ions  is  then  |H*rpondictilnr  to  the  direct  ion  of  propagation.)  h'or  this 
case,  it  is  convenient  to  return  to  equation  .’?(>  so  that  with  \  0,  we  have 


j  A  4> 

ros  i/o.sin; 

>  </  ;> 

which  for  A  9  •:  1  gives 

|F(A<p)|  =  S1.'}  Kv  A  9  (40) 

ky  A  9 

For  A  9  •  '  1 .0,  A  9=  sin  Ac  so  AA  9  =  ky=  2  tt  Av  (see  fig.  9).  There¬ 
fore  the  v  coherence  provides  a  measure  of  the  dimension  perpendicular 
to  the  direction  of  propagation,  i.e.,  the  "long-crestedness”  of  the  waves. 

The  v  coherence  for  separation  in  the  v  direction  is  shown  plotted 
as  the  white  curves  in  figure  10  and  the  right-hand  scale  should  be  used. 

It  is  immediately  seen  that  the  coherence  falls  off  much  more  rapidly  (for 
a  given  beam  width)  lor  station  separation  perpendicular  to  the  direction 
of  propagation  than  for  separation  along  the  propagation  direction. 

The  general  situation  when  E(co)  is  independent  of  9  and  A  9-  1.0 

is  summarized  in  figure  11  where  the  \  conerenees  for  arbitrary  diree-  - 
tions  of  separation  are  shown  for  an  assumed  station  separation  of 
L  A  .  V?. 

The  other  limiting  case  that  is  interesting  to  consider  is  that  of 
nondirectional  radiation  (uniform  in  all  directions)  so  that  A 9=  tt.  Then 
we  see  from  equation  -‘JO  that 

|fU'v>|  —  f  cos  (/.’v sin 9)  c/9  “  jj0(/'V)|  (41) 

TT  CO  1 

The  v' coherence  then  falls  off  with  station  separation  ( v  L>  as  the  '■■Jn(hy)\ 
curve  shown  in  figure  12.  The  coherence  becomes  negligible  in  all  direc¬ 
tions  will'll  the  separation  is  about.  0.4  wavelength. 

if  the  recorders  are  located  very  near  a  large  source  region  it  is 
of  interest  to  consider  the  integral  taken  over  the  whole  half-plane  that 
includes  the  source  t— tt  2  9  tt  2 > .  Then 


1  f-n  ~ 

l‘  /c yf  ■  —  /  cos  ( Av  sin  exp  (  —  ik\  cos  9)  1/9  (4 

TT  J  r T  2 

h  or  sopara!  1011  aiong  the  \  direction  <  \  0*  the  integral  is  again  -A, 'Ay* 

since  on  ehanging  the  \aiiablo  9  to  where  < u  2»-q> 


1 


I 


Figure  12.  Dependence  ol'\  coherence  on  recorder  separation  and 
wavelength  when  the  integration  extends  over  all  angles  (isotropic 
ease).  \  coherence  then  decreases  equally  for  all  directions  of 
separation  according  to  J0  (/■•>).  When  integration  is  over  Lhe  whole 
half-plane  of  the  source,  \  coh  lulls  off  according  to  |F(A.t)| 

for  separation  in  the  x  direction  and  according  to  ;F Uty)!  for 
separation  in  the  y  direction. 


FU fv>  ~  -  cos  (/I’vcosy)  </ v  -  f/0(/dv) 

TT  /()  ' 


For  separation  along  the  x  direction  (y  =  0) 


Fl/c.v)  =  I 


cTT/2  .  /-TT.  2 

I  /  ^  cos  ( /e v  cos  <p)  <7 9  -  —  /  ^  s  i  n  ( A’ a  cos  9^  (/  9  ( 


By  a  transformation  similar  to  that  above  (i.e.,  letting  v  -  (tt  2) -9  the 
first  integral  is 


1  rn 

—  /  cos  (A.vsin  v  >  (/y  -  J„(A.\) 


I 


The  second  integral  is  more  complicated.  Again  transforming  to  a  variable 
ranging  from  limits  ofO  to  tt, 


i  r 2 

TT  J-T T/i 


sin  (kx cos 


sin  (kx  sin  y )  cly 


The  integrand  sin  (kx  sin  y)  can  be  represented  as  a  series  (see,  for  ex¬ 
ample,  ref.  13)  so  that 


_  2 
y  — 

TT 

2 

Ji  n 4  . (  kx) 

X  sin 

n=  0 

_  2 

OG 

Jin  +  i  (  kx)  | 

'  2  \1 

TT 

z 

n  -0 

An*  1/ 

II 

=t  I4- 

J,  (kx) 
o 

(46) 


So  for  separations  in  the  x  direction  F(kx )  is  simply  evaluated  from 


F(kx)  -  J0  ( kx)  -  —  i 


jam  . 

3  o 


For  kx  '  .  2  tt  only  a  few  terms  are  needed  and  the  Bessel  functions  are 
conveniently  tabulated  in  many  places.  F(kx) j  is  shown  in  figure  12  for 
kx  ■  2  tt. 


Sample  Problem  Applying  Cross-Spectrum  Analysis  to 
Spatially  Separated  Recordings 

Figure  13  shows  the  configuration  of  an  experiment  reported  by 
Gossard  and  Paulson  14.  Radio  waves  at  a  frequency  of  44.3  kHz  were 
reflected  off  the  lower  ionosphere  and  received  at  a  triangle  of  spaced 
receivers.  The  receiver  sites  are  indicated  by  MW,  OR,  and  GL  and  their 
separations  in  kilometers  are  shown  on  the  figure.  It  is  assumed  that  the 


39 


* 


St 

I 

m 

■m 

'S 

=1 

= 

4 


movement  of  ionospheric  irregularities  in  electron  density  causes  a  simi¬ 
lar  movement  of  the  diffraction  pattern  of  reflected  radio  waves  received 
at  the  surface  of  the  earth.  The  observed  movement  of  the  diffraction 
pattern  is  then  compared  with  ionospheric  wind  measurements  obtained  by 
photographing  the  drift  of  chemiluminescent  trails  produced  by  contami¬ 
nants  injected  into  the  lower  ionosphere  by  a  gun-launched  projectile15 
located  near  midpath. 
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Figure  18.  Perspective  drawing  showing  geometry  of  multiple  receiver 
experi iiu  ni  to  iiH’iisuri'  movement  of  irregularit  n*s  of  electron  density  m 
the  ionosphere •  Length  units  in  MW.  OK.  (II.  receiver  triangle  arc  kilo- 
ineters.  Transmitter  is  near  Sent  niel .  Arizona.  Lata  are  analyzed  for 
ionospheric  nave  motions.  using  cross-spectrum  technii|ue. 
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Sample  records  of  amplitude  fluctuations  observed  at  the  three 
receiver  sites  on  13  and  14  June  1966  are  shown  in  figure  14,  together 
with  height  contours  of  the  500-mb  pressure  surface  two  days  earlier.  The 
time  lag  is  roughly  in  accord  with  the  theoretical  time  required  for  energy 
to  flow  from  the  troposphere  to  the  ionosphere  for  typical  atmospheric 
wind  and  temperature  distributions.  Therefore  it  is  of  interest  to  inquire 
whether  the  direction  of  propagation  of  the  diffraction  pattern  and  its 
'‘beam  width”  are  compatible  with  such  a  source. 

A  cross-spectrum  analysis  of  the  records  of  the  three  stations 
was  performed,  and  it  was  found  that  the  velocity  of  propagation  of  a 
prominent  spectral  line  with  a  period  of  5.8  minutes  was  47  mps  toward 
a  direction  14  degrees  west  of  north.  The  wave  crests  were  therefore  only 
about  12  c  from  alignment  with  the  MW-GL  line  of  separation  of  receivers. 
The  indicated  coherence  between  MW  and  GL  for  this  spectral  line  was 
about  0.44.  The  wavelength  of  this  wave  component  was  17  km  and  the 
station  separation  was  21  km,  so  L  ' 1.24.  Figure  10  shows  that  A  <p 
is  therefore  about  0.3  radian  so  that  the  transverse  wavelength  is  about 
triple  that  normal  to  the  crests.  We  therefore  conclude  that  the  direction 
of  propagation  of  the  5.8-minute  spectral  line  is  compatible  with  the  ob¬ 
served  position  of  the  tropospheric  disturbance  and  the  '’beam  width” 
indicated  by  the  coherence  analysis  is  appropriate  to  an  extended  source. 


I 


Figure  14.  Sample  reeord  of  radio  amplitude  fluctuations  (top)  produced 
!>.v  wave  motions  in  the  lower  ionosphere.  Radio  frequency  is  44.3  kHz, 
time  is  MST.  and  amplitudes  are  normalized  to  the  record  average.  The 
"outlier  map  (  bottom)  shows  height  contours  of  the  500-nib  pressure 
surface  tvo,  days  earlier.  The  white  triangle  shows  the  location  of  the 
expi  muelilal  site. 
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APPENDIX  A:  CONVENIENT  NUMERICAL  FILTERS 

One  of  the  simplest  mathematical  filters  is  a  running  average. 
Consider  the  effect  of  a  running  average  on  the  Fourier  component 

.  2  TT?  .  . 

A  =  A  cos  ——  ( A- 

When  the  running  average  extends  from  -  ta/ 2  to  +  tQ  /2  the  average  is 
taken  symmetrically  about  the  maximum  of  a  cosine  curve,  and  we  have 
for  the  running  average,  a 


-  .4  f‘a  2  2  TTt 

A  -  —  /  COS  - 

laJ-W  2  T 


(A-2) 


Integration  of  (A-2)  gives  the  filter  response  F<T)  as 


FIT)  „  -1  sin  ^ 


(A-3> 


The  power  at  each  period  T  is  reduced  according  to 


FMT)  •= 


( A—4) 


Figure  A1  shows  FJ(7’)  or  the  filter  response  of  the  simple  running  average 
as  the  white  curve.  The  high  frequencies  are  suppressed  but  FJ(T)  contains 
large  side  lobes,  the  first  of  which  has  an  amplitude  of  about  one-tenth  of 
the  main  lobe.  This  filter  is  undesirable  because  side-lobe  frequencies 
would  be  accented  in  a  spectral  analysis. 

A  more  desirable  filter  results  from  a  triangularly  weighted  average 

given  by 


4V'2(l-t)- T“" 


(.4-5) 


1 


with  a  filter  response 


F(T)  = 


•(f)'- 


(A-6) 


F 2  (T)  from  equation  A-6  is  shown  as  the  black  curve  in  figure  Al.  The 
function  has  smaller  side  lobes  than  that  of  the  filter  given  by  equation 
A-4,  but  the  frequency  cutoff  is  less  abrupt. 

The  low-pass  filter  (eq.  A-6)  can  be  used  as  a  high-pass  filter  by 
subtracting  the  filtered  data  from  the  unfiltered  (raw)  data.  The  response 
of  this  high-pass  filter  is  given  in  figure  A2  calculated  using  A-7. 

The  filter  can  be  written  for  digitized  data  as  the  right-hand  term  of  the 


following  equation  where  y,  =  x,  -  x  and  m 


2\A  t  / 


yi=Xi~.\Xi-m  +  l+2xi.m  +  2  +  +  ~  1>  X,  .  x  +  mxt  +  (m  -  1)  X .  +  j 


(A-7) 


+  •  •  • 


+  2xI  +  m.2  +  xi  +  m.1 


where  i  is  the  i-th  value  in  the  time  series  and  m  is  the  half-width  of  the 
triangularly  weighted  averaging  interval.  The  sample  size  is  reduced  by 
2m  (m  points  are  lost  both  at  the  beginning  and  the  end  of  the  sample)  and 
the  mean  of  the  filtered  data  is  zero.  Other  filters  can  be  designed  but 
those  given  by  equation  A-4  or  A-6  are  convenient  and  adequate  for 
many  purposes. 

The  data  shown  in  figure  2  were  filtered  using  an  m  of  100  to 
remove  the  lower  frequencies.  This  filter  reduced  the  amplitude  of  the 
annual  cycle  of  fE  to  about  0.05  of  its  original  value  while  leaving  com¬ 
ponents  having  periods  less  than  100  days  essentially  unchanged.  The 
filtered  data  are  shown  as  the  bottom  record  of  figure  2.  Because  the 
first  and  last  100  data  values  were  eliminated  by  the  filter,  the  raw  data 
extended  100  days  before  and  100  days  after  days  numbered  0  and 
180,  respectively. 

Examination  of  the  filtered  data  in  figure  2  clearly  shows  that  the 
higher  frequencies  are  almost  unaffected  by  the  filter.  The  reduction  in 
the  amplitude  of  the  component  having  the  period  of  the  annual  cycle  is 


A-3 
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obvious  but  no  obvious  reduction  is  noted  in  periods  of  100  to  200  days. 

Analysis  of  the  filtered  fE  should  not  extend  to  periods  much 
greater  than  100  days.  This  study  used  a  sample  size  of  2722  (days). 

For  the  Fourier  component  having  a  period  of  100  days,  /'£  has  27  possible 
cycles  in  the  record.  Twenty-seven  cycles  are  considered  adequate  and 
are  considered  to  be  within  the  large  sample  category. 
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APPENDIX  B:  USEFUL  FORMULAE  FOR  DI6ITAL 
COMPUTATION  OF  SPECTRA 


A  u  to-C  o  v  a  r  i  a  nee : 


i  -•  \  -  p 


Ru  <p> 


Y_p  ^  (.v,  -  x, ),-  < Xj  - x2 ) 

<  - 1 


where 


i  :\'  -  P 


/  .=  V 


l  V  l  V 

*1  =  - -  >  A,-,  -Vj  - -  >  Xj, 

iN”P  A-p 


<  - 1  -  p 


Cross-Covariance  Function: 


/  A  -  p 


Ru<p) 


i  .  i 


where 


a  - 
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i  \  -  p 


/  =  \ 


A  -  P 


i  -  -a  i 


^  f  -*  p 


/  -  1  i  •  4  )  p 

Power  Spectrum  ( per  cycle  per  unit  time): 


p  .  m 


Elf)  -4 A  fS  ^  /?j,(p)  cos  where  5 


w 


p  n 


1  2  for  p  -  0 ,  m 
1  for  0  p  ■'  m 


Power  Spectrum  (Includes  Hanning  filter  window)! per  cycle  per  unit  time): 
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Cross-Power  Spectra: 


C(f)  .  2A/ 
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Bandwidth  of  Analysis: 


A/  -  (2mA  /)"' 

Highest  Frequency  in  Digital  Analysis  (Nyguist  frequency): 

-Coherence  (/)  ^  L2,n  +  Q  (Q 
E, </')  E2(/) 

Phase  Cross  Spectrum: 


tan  cot  =  tan  2  tt/‘t  -  21l! 

C  <  /  > 
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This  document  attempts  to  bring  together  a  description  of  the  methodology  of  spectrum 
and  cross  spectrum  analysis  which  is  scattered  in  the  scientific  literature.  The 
digital  formulation  required  by  most  engineers  and  applied  physicists  are  given. 
Elementary  derivations  of  the  basic  formulae  are  given.  The  various  errors  that  may 
contaminate  the  analysis  are  discussed  together  with  the  confidence  limits  applicable 
to  samples  of  finite  length.  Examples  ar'1  given  illustrating  the  application  and  the 
step-by-step  procedure  for  obtaining  spectra  and  cross  spectra.  Two  very  different 
applications  of  cross-spectrum  analysis  are  described,  and  the  significance  of  cross¬ 
power  and  coherence  is  discussed. 
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